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1374-SPL C 

METHOD AND APPARATUS FOR 
SEGMENTING SMALL STRUCTURES 
IN IMAGES 

Cross-Reference to Related Applications 

The present indention is related to the subject matter of provisional patent 
application Serial No. 60/084,125 filed on May 4, 1998, the entire disclosure of 
which is incorporated herein by reference. 

Field of the Invention 

The present invention relates to data processing of intensity data arranged 
in a multidimensional array. More particularly, the invention relates to a method, 
an apparatus, and computer program products for rapidly segmenting 
multidimensional intensity data by which points in one or more small structures 
contained in the data are labeled. 

Background of the Invention 

Digital imagery and other multidimensional digital arrays of intensity are 
routinely collected using digital sensors and arrays of charge coupled devices 
(CCDs). The resulting data arrays are analyzed to determine patterns and detect 
features in the data. For example, color images of a battle scene are analyzed to 
detect targets, and radiographs and sonograms of human and animal bodies are 
examined to detect tumors and other indications of injury or disease. As the 
number and complexity of these digital data arrays to be analyzed increase or the 
time required to perform the analyses decreases, automated and machine assisted 
analysis becomes more critical. Some statistically based automated procedures 
for detecting features in a multidimensional array are adequate when the feature 
encompasses many points in the array, i.e. when the feature is large, but fail to 
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perform well as the feature to be detected becomes small. Some procedures 
perform well when tuned to a particular problem through experimental adjustment 
of many parameters, but such tuning may place an undue burden on time limited 
or experience limited personnel. Typical problems encountered with such 
5 automated analysis of small structures in multidimensional arrays are illustrated 
for the case of automatic detection of microcalcification candidates in 
mammograms. 

Breast cancer has the highest incidence among all cancer types in 
American women, causing 1 woman in 8 to develop the disease in her lifetime. 

10 Every year, about 182,000 new cases are diagnosed with breast cancer and about 
46,000 women die of this disease. The 5-year survival for women with breast 
cancer improves significantly with early diagnosis and treatment. To enable early 
detection, the American Cancer Society (ACS) recommends a baseline 
mammogram for all women by the age of 40, a mammogram approximately every 

15 other year between the ages of 40 and 50, and a mammogram every year after the 
age of 50. It is possible that the volume of mammography will become one of the 
highest among clinical X-ray procedures since more than 30 million women in the 
US are above the age of 50 and 41% are known to follow the ACS guidelines. 

Besides the volume problem, an additional difficulty of early detection of 

20 breast cancer in mammograms is the subtlety of the early signal. A 
microcalcification cluster, an early sign of breast cancer that may warrant biopsy, 
is commonly defined as three or more microcalcifications present in 1 cm 2 on a 
mammogram. These clusters are often difficult to detect due to their small size 
and their similarity to other tissue structures. The width of an individual 

25 microcalcification is less than 2 mm. The etiology of microcalcifications includes 
lobular, ductal or epithelial hyperplasia, secretion of calcium salts by epithelial 
cells, adenosis, as well as calcification of necrotic debris due to carcinoma. Up to 
50% of breast cancer cases exhibit microcalcification clusters, and 20-35% of 
clusters in the absence of a mass are related to malignant growth. In many cases a 

30 cluster is the first and only sign that allows timely intervention. 



The increasing pressure to interpret large numbers of mammograms and 
the subtlety of many early signs increase the likelihood of missing breast cancer. 
A reliable automated system that indicates suspicious structures in mammograms 
can allow the radiologist to focus rapidly on the relevant parts of the mammogram 
and it can increase the effectiveness and efficiency of radiology clinics. In the 
detection of breast cancer, false negatives may cause a delay in the diagnosis and 
treatment of the disease while false positives cause unwarranted biopsy 
examinations. Therefore, both sensitivity and specificity need to be maximized, 
with a relatively higher priority on sensitivity, which has a more vital role. 

A common approach used for detecting microcalcifications in 
mammograms starts by segmenting candidate structures and subsequently 
applying feature extraction and pattern recognition to distinguish 
microcalcifications from background tissue among the candidates. In this process, 
segmentation plays an essential role since the quantitative features that represent 
each candidate structure, such as size, contrast, and sharpness, depend on the 
region indicated by segmentation. Furthermore, to process all possible candidate 
structures, a considerably large number of background structures need to be 
segmented, making fast segmentation desirable. 

Several techniques for segmentation have been applied to 
microcalcifications. One segmentation technique is based on local thresholding 
for individual pixels using the mean pixel value and root mean square (rms) noise 
fluctuation in a selected region around the thresholded pixel. The threshold for a 
pixel is set as the mean value plus the rms noise value multiplied by a selected 
coefficient. A structure is segmented by connecting pixels that exceed the 
threshold. Both parameters that have to be selected, size of region and threshold 
coefficient, are critical to this method. If a microcalcification is close to another 
microcalcification or bright structure, the window used to compute the rms noise 
value around the first microcalcification will include the other bright structures, 
and the noise rms may be overestimated, thus setting the threshold too high. On 
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the other hand, if the selected region is too small, it will not contain sufficient 
background pixels when placed on large microcalcifications. 

Such a window size needs to be selected in a second segmentation 
algorithm as well, where local thresholding is used by setting a threshold for small 
5 square sub images. The threshold is based on an expected bimodal intensity 
distribution in a window of selected size that contains the sub-image to be 
segmented. If the distribution is not bimodal, then the threshold is set by using 5 
different positions of the window each containing the sub-image to be segmented. 
The existence of a bimodal distribution in at least one window is essential for this 
10 algorithm. 

Other segmentation methods start with seed pixels and grow a region by 
adding pixels. They also require selection of a window size and threshold 
parameters. The localized implementation of region growing depends on the 
selected window size and the threshold for absolute difference in gray level 
15 between the seed pixel and a pixel to be added to the region. 

One segmentation algorithm uses several steps that include high-pass 
filtering, difference of Gaussinan filtering, four computations of the standard 
deviation of the image, a smoothing, an opening, as well as an iterative thickening 
process with two erosions, two intersections and a union operation in each 
20 iteration. More than ten parameters have to be selected, including widths of 
Gaussian distributions, threshold coefficients, and diameters of morphological 
filtering elements. 

A segmentation algorithm that operates without parametric distribution 
models, local statistics windows, or manually adjustable thresholds is desirable. 
25 A segmentation method that is fast is also important. Up to 400 films per 

day are routinely screened in busy radiology clinics. The automated analysis does 
not have to be applied on-line; however, it may be difficult to process large 
numbers of mammograms overnight if algorithms are not fast enough. Because 
the segmentation algorithm has to segment all candidate structures that may 
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potentially be microcalcifications, its speed is especially relevant. Each film may 
have several thousand candidate structures that must be segmented. 

The multi-tolerance segmentation algorithm of Shen et al. (L. Shen, et al. 
"Detection and Classifications of Mammpgraphic Calcifications," International 
5 Journal of Pattern Recognition and Artificial Intelligence, vol. 7, pp. 1403-1416, 
1993), does not use statistical models for local statistics, and its threshold is set 
automatically. This multi-tolerance, region growing approach uses a growth 
tolerance parameter that changes in a small range with a step size that depends on 
the seed pixel. The structure of interest is segmented multiple times with varying 

10 tolerance parameters, and in each segmentation, a set of three features is 
computed. The normalized vector differences in the feature set between 
successive segmentations are calculated and the segmentation with minimal 
difference is selected as the final one. 

The active contours model of Kass et al. (Kass, M. et al. "Snakes: Active 

15 Contour Models," International Journal on Computer Vision, pp. 321-331, 1988), 
also provides segmentation without parametric statistical data models or windows 
for local statistics, but does rely on several user selected parameters that place 
some burden on the user. It has been used successfully to determine the 
boundaries of tissue structures in data such as ultrasound and MRI images of the 

20 heart, and MRI images of the brain, but it has not been applied to the 
segmentation of microcalcifications. The active contours model starts with an 
initial contour placed near the expected boundary and moves the contour 
iteratively toward the boundary by minimizing an energy function. The contour is 
modeled as a physical flexible object with elasticity and rigidity properties. Its 

25 dynamics, dictated by the balance between these internal properties and external 
forces that depend on the image data, satisfy . the Euler equations and minimize the 
corresponding energy function. An active contour that is initiated as a closed 
curve remains so during iterations and its smoothness can be adjusted by the 
choice of parameters. 
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What is needed is a segmentation method and apparatus without statistical 
models, local statistics, or thresholds to be selected manually, and with 
significantly lower computational complexity compared to the multi-tolerance and 
active contours methods, for enhanced speed. 
5 In particular, what is needed is a method and apparatus to segment pixels 

in an image, such as a mammogram, containing a plurality of extra dark or extra 
bright objects just a few pixels in extent, that gives edges similar to those selected 
by an expert, but does so with fewer computations and with fewer manually 
adjustable parameters than conventional segmentation methods and equipment. 

10 Summary of the Invention 

Therefore it is an object of the present invention to provide segmentation 
for small features in multidimensional data which defines small feature edges that 
correspond closely to those selected by an analyst but does so with less 
complexity than the above known methods. 
15 It is another object of the present invention to provide a data processing 

apparatus that more rapidly provides small feature edges that correspond closely 
to those selected by an analyst. 

It is another object of the present invention to provide computer program 
products that more rapidly provide small feature edges that correspond closely to 
20 those selected by an analyst. 

It is another object of the invention to identify microcalcifications in a 
mammogram. 

These and other objects and advantages of the present invention are 
provided by a method for segmenting a small feature in a multidimensional digital 
25 array of intensity values in a data processor. Each small feature includes a local 
intensity extreme, such as an intensity maximum. An edge metric is computed 
along each ray of a plurality of multidimensional rays originating at the local 
intensity extreme. A multidimensional edge point is identified corresponding to a 
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maximum edge metric on each ray. Every point on each ray from the local 
extreme to the ray edge point is labeled as part of the small feature. The labeling 
is then spread to an unlabeled point following a hill climbing procedure requiring 
that the unlabeled point be adjacent to a labeled point, have a similar or more 
5 extreme intensity than the labeled point, and be closer than the labeled point to the 
local extreme. 

In another embodiment, the multidimensional array is a digital image, and 
each point is a pixel. In another embodiment, the digital image is a digitized 
mammogram and the small feature is a microcalcification candidate. In the latter 

10 embodiment, microcalcification candidates are satisfactory segmented in fewer 
operations than with conventional segmentation methods. 

In another aspect of the invention, a data processing apparatus segments a 
small feature in a multidimensional digital array of intensity values. The 
apparatus includes an input for inputting a plurality of intensity values arranged 

15 along regular increments in each of a plurality of dimensions and a memory 
medium for storing the plurality of intensity values as a multidimensional digital 
array. The apparatus includes a processor configured to detect a local intensity 
extreme in the multidimensional digital array, to identify points along a plurality 
of rays originating at the local intensity extreme, and to identify one ray edge 

20 point on each ray. The ray edge point is associated with a maximum edge metric 
along the ray. The processor is also configured to label the points in the array that 
are part of the small features. Each point on each ray from the local intensity 
extreme to the edge point is labeled, as is an unlabeled point adjacent to a labeled 
point if the unlabeled point has a more extreme intensity than the labeled point 

25 and the unlabeled point is closer than the labeled point to the local extreme. 
Labeling continues until no more unlabeled points can be labeled. The apparatus 
also includes an output for providing the labeled points for subsequent processing. 

In another aspect of the invention, a computer program product is provided 
for segmenting a small feature in a multidimensional array of intensities using a 

30 computer. The computer program product includes computer controlling 
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instructions for configuring a computer to compute an edge metric along each ray 
of a plurality of multidimensional rays originating at a local intensity extreme. 
The instructions also identify a ray edge multidimensional point corresponding to 
a maximum edge metric on each ray. The program also labels every point on each 
5 ray from the local extreme to the ray edge point, and then labels an unlabeled 
point if the unlabeled point is adjacent to a labeled point and the unlabeled point 
has a more extreme intensity than the labeled point, and the unlabeled point is 
closer than the labeled point to the local extreme. In one embodiment, the 
instructions are stored in a computer readable memory device. In another 
10 embodiment, the instructions are transmitted as electronic signals on a 
communications line. 

The foregoing and other features, aspects and advantages of the present 
invention will become more apparent from the following detailed description of 
the present invention when taken in conjunction with the accompanying drawings. 

15 Brief Description of the Drawings 

The preferred and example embodiments of the present invention are 
described with reference to the Drawings in which: 

Figure 1A is a perspective view of the external features of a computer 
apparatus suitable for one embodiment of the present invention. 
20 Figure IB is a block diagram of a computer apparatus that can be 

configured according to one embodiment of the present invention. 

Figure 1C is a perspective view of a sample memory medium for storing 
instructions to configure a computer according to another embodiment of the 
present invention. 

25 Figure ID is a block diagram of a network that can transmit electronic 

signals that configure a computer according to still another embodiment of the 
present invention. 
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Figure 2A is a flow diagram for a method according to yet another 
embodiment of the present invention. 

Figure 2B is a flow diagram following step 270 of Figure 2A according to 
a further embodiment of the present invention. 
5 Figure 2C is a flow diagram for details of step 260 of Figure 2 A according 

to still another embodiment of the present invention. 

Figure 2D is a flow diagram for an alternative detail for step 260 of Figure 
2A according to yet another embodiment of the present invention. 

Figure 3 is a schematic diagram of a local maximum, rays and edges that 
1 0 results from steps 2 1 0 through 250 of Figure 2. 

Figure 4 is a schematic diagram of a local maximum, a labeled pixel, 
adjacent pixels, and a reference line according to one criteria for one embodiment 
of step 260 of Figure 2. 

Figure 5 is a schematic diagram of a local maximum, a labeled pixel, and 
15 an adjacent pixel according to a criteria for another embodiment of step 260 of 
Figure 2. 

Figures 6A-6D are gray scale photographs showing an actual intensity 
maximum as originally provided and then superposed with labeled pixels after 
three stages of the method of Figure 2 according to the present invention. 
20 Figures 7A-7D are gray scale photographs showing three actual intensity 

maxima as originally provided and then superposed with labeled edge pixels after 
segmentation based on two conventional methods and the preferred embodiment 
of the present invention. 

Detailed Description of the Preferred Embodiment 

25 The principles of the present invention will be described next, detailed in 

terms of preferred and example embodiments with reference to the accompanying 
drawings. Whenever possible, the same reference numbers will be used 
throughout the drawings to refer to the same or like parts. 
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The explanations of the detailed embodiments are by way of example only 
and are not meant to limit the scope of the invention. The invention applies to 
identifying small structures in any multidimensional array of regularly spaced 
intensity values. Here intensity is used in a generic sense representative of 
5 measured data values in general, and is not confined to density of optical energy. 
Examples of such multidimensional arrays include gray-scale digital images in 
which intensity values are regularly spaced in two dimensions, often called rows 
and columns or y and x, such as the mammogram described in the preferred 
embodiment. In this kind of arrangement, each digital image element is a picture 
10 element called a pixel. Elevation maps are two dimensional arrays of height data, 
where height is the "intensity." Other examples of multidimensional arrays 
include color images which can be represented as three-dimensional arrays of 
intensity where the third dimension is color. Typically, the array would have 
intensity at only three points in the color dimension, for example, a red intensity, a 
15 blue intensity and a green intensity. Gray-scale video clips can also be considered 
three-dimensional arrays, where each video image frame is two-dimensional and 
the third dimension is time. By the same token, color video clips can be 
considered four-dimensional where the four dimensions are row, column, color 
and time. Other examples include medical imagery where two-dimensional cross 
sections of a human body are assembled at several positions from head to toe. In 
this case the third dimension is height through the subject. By extension, such 
three-dimensional looks can be repeated at uniform intervals of time, making time 
the fourth dimension. Thus the descriptions that follow apply not only to gray 
scale images of the preferred embodiment, but to multidimensional arrays of 
digital data. 

A multidimensional point in a multidimensional digital array is located by 
the index of the point in each of the dimensions. Letting D represent the number 
of dimensions, the location of a multidimensional point P in a multidimensional 
array can be specified uniquely by a set containing D indexes as coordinates, {I„ 
h> Ij. ■ • -Id)- Where there are only two dimensions, it is common to refer to l x as 
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the x coordinate and to refer to I 2 as the y coordinate. There is an implied limit to 
the number of allowed positions in each dimension of a finite array. Letting L { 
represent the maximum number of locations in the i-th dimension of the digital 
data array, each index can vary from one to Lj, inclusive. That is: 
5 l<Ii<L, (1) 

The distance, d, between any two multidimensional points, P a and P b , with 
different indices {a l5 a^ . . .a D } and {b x , b 2 , b 3 , . . .b D }, can be defined as the 
square root of the sum of the squares of the differences in their indices. That is, 
d(P a , P b ) = d(P( ai , a,, . . ., a D ), P(b„ b 2 , . . ., b D )) 
10 = V [OvaJ 2 +(b 2 -a 2 ) 2 + . . . + (b D -a D ) 2 ] (2) 

The intensity, f, varies with position in the multidimensional array and 
may be represented by the symbol f(P). The intensity f at each multidimensional 
point can be a single value, also called a scalar quantity. Alternatively, the 

15 intensity can be a vector of several values, e.g., f(P) = {fl(P), f2(P), f3(P)}. For 
example, the three-color image can be treated as a three-dimensional array or can 
be treated as a two dimensional image with a three element vector intensity. In 
this terminology, the vector elements of the intensity are not used in the 
calculation of distance using Equation 2. Instead, the magnitude of intensity at 

20 point P could be any vector magnitude convention such as the square root of the 
sum of the squares of the vector components or the sum of the absolute values of 
the vector components. Similarly, the difference in intensity between two points 
P a and P b would be given by the magnitude of the difference in the components 
using any conventional method. 

25 Thus, though the preferred embodiment is described in which the digital 

data array is an image having two dimensional pixels, each pixel having a scalar 
image intensity, the method can readily be extended to multiple dimensions using 
the above relationships. In the following, each pixel P has a first coordinate 
represented by x and a second coordinate represented by y and an intensity 



12 

represented by f(P) or f(x,y). Separate pixels are designated by separate 
subscripts. 

Though the invention applies to any imagery, the preferred embodiments 
segment two-dimensional images with a gray-scale intensity representative of a 
5 mammogram. Other two dimensional imagery which the present invention can 
segment include imagery of military scenes in which the intensity is responsive to 
the presence of targets of a firing system, such as vehicles to be fired upon by a 
missile. 

The invention is related to finding small objects in a multidimensional 
10 array. In this context small means objects affecting the intensity in several points 
in one dimension of the array but not many thousands of points in each 
dimension. Other, statistical and textural segmentation procedures, for example, 
are expected to be more useful as the number of points in a feature increases. It is 
characteristic of microcalcifications in mammograms and distant targets in 
15 military scenarios that only several pixels are contained in the object to be 
segmented. It is also anticipated that many other features to be detected in 
radiographs and sonograms of biological bodies also involve only several pixels. 
The present invention is expected to perform especially well for these 
applications. 

20 The methods and procedures discussed herein are intended to be 

performed by data processing systems or other machines. Though described in 
terms that can be interpreted to be performed by a human operator, such 
performance is neither required nor likely to be desirable. Multiple tedious 
computations with high accuracy are required that are unsuitable for practical 

25 implementation by human beings. Also, the invention can be implemented in 
computer or other hardware, the structure of which is evident from the following 
descriptions. 

Also herein, the procedures will be described as the manipulation of 
values, symbols, characters, numbers, or other such terms. Though such terms 
30 can refer to mental abstractions, herein they are used as convenient expressions 
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for physical signals such as controllable chemical, biological, and electronic and 
other physical states that can be used to represent the values, symbols, characters, 
numbers, or other such terms. 

Figure 1A illustrates a computer of a type suitable for carrying out the 
5 invention. Viewed externally in Figure 1A, a computer system has a central 
processing unit 100 having disk drives 11 OA and HOB. Disk drive indications 
11 OA and HOB are merely symbolic of a number of disk drives that might be 
accommodated by the computer system. Typically these would include a floppy 
disk drive such as 1 10A, a hard disk drive (not shown externally) and a CD-ROM 

10 drive indicated by slot 1 10B. The number and type of drives vary, typically, with 
different computer configurations. The computer has a display 120 upon which 
information is displayed. A keyboard 130 and mouse 140 are typically also 
available as input devices. 

Figure IB illustrates a block diagram of the internal hardware of the 

15 computer of Figure 1A. A bus 150 serves as the main information highway 
interconnecting the other components to the computer. CPU 155 is the central 
processing unit of the system, performing calculations and logic operations 
required to execute programs. Read-Only-Memory 160 and Random- Access- 
Memory 165 constitute the main memory of the computer. Disk controller 170 

20 interfaces one or more disk drives to the system bus 150. These disk drives may 
be floppy disks drives, such as 173, internal or external hard drives, such as 172, 
or CD-ROM or DVD (digital video disk) drives such as 171. A display interface 
125 interfaces a display 120 and permits information from the bus to be viewed 
on the display 120. Communications with external devices can occur over 

25 communications port 175. 

Figure 1C illustrates an exemplary memory medium which can be used 
with drives such as 173 in Figure IB or 11 OA in Figure 1A. Typically, memory 
media such as a floppy disk, or CD-ROM, or DVD, will contain the program 
information for controlling the computer to enable the computer to perform its 

30 functions in accordance with the invention. 
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Figure ID is a block diagram of a network architecture suitable for 
carrying data and programs over communication lines in accordance with some 
aspects of the inventions. A network 190 serves to connect a user computer or 
client computer 1 10 with one or more servers such as server 195 for the download 
5 of program and data information. A second user on a second client computer 100' 
can also connect to the network via a network service provider, such as ISP 180. 

In general, small objects in images may have an intensity level that is 
either lower or higher than a surrounding background. An intensity maximum 
with levels higher than the background is called a local maximum, and an 

10 intensity minimum with intensity levels below the background is called a local 
minimum. Both maximum and minimum are encompassed by the term intensity 
extreme. Thus, in general, the target objects in an image or multi-dimensional 
array encompass intensity extremes. Both are capable of being segmented 
according to the present invention. For the sake of serving as an example, the 

15 following description generally considers the preferred embodiment in which 
microcalcifications are evident as local maxima in intensity, and the method will 
be called a hill climbing method; however, segmenting a local minimum is also 
anticipated using the hill climbing method. In the following discussion, when a 
first point has an intensity equaling the intensity of the local extreme or between 

20 the intensity of the local extreme and the intensity of a second point, the first point 
is said to have a more extreme intensity than the second point. 

Figure 2A shows the method according to one embodiment of the present 
invention. A local brightness maximum, characteristic of a microcalcification, is 
identified at pixel P 0 in an image at step 210. Next, a plurality of rays is defined 

25 that emanate from that local maximum pixel P 0 as illustrated in step 220. Figure 3 
illustrates five sample rays 320 emanating from a local maximum 310. Referring 
again to Figure 2A, an edge metric is computed for each pixel along each ray in 
step 230. Then in step 240, a ray edge pixel on the ray is identified based on a 
maximum edge metric. Then the pixels on the ray from the local maximum to the 

30 ray edge pixel, inclusive, are labeled as belonging to the object or feature in step 
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250. Additional pixels belonging to the feature are labeled if they are adjacent to a 
labeled pixel and if the unlabeled pixel satisfies intensity and distance criteria 
described later. These criteria implement the unique hill climbing procedure of the 
present invention. This growth of labeled pixels is indicated by step 216. In step 
5 270, every unlabeled pixel next to a labeled point is examined using the criterion 
in step 260 until no further points can be labeled. 

Figure 2B shows steps that follow step 270 in another embodiment of the 
present invention. Here each of the labeled pixels is checked in step 275 and those 
labeled pixels adjacent to an unlabeled pixel are relabeled as an edge pixel of the 

10 small feature. This completes the labeling associated with one of the small 
features in the image; and, in step 280, control is returned to step 210 until no 
local maximum remains unlabeled or unsegmented in the image. In yet another 
embodiment of the invention, small features identified in the image can be joined 
in step 285 if those pixels are within a joint distance. Additional detail regarding 

15 the steps shown in Figures 2 A and 2B are provided with reference to Figures 2C 
through 5. 

According to the present invention, the segmentation is based on the 
experience that, in a given array, the edge of a small feature to be segmented is a 
closed contour around a local intensity extreme pixel P 0 . In the preferred 

20 embodiment, the local intensity extreme is selected as the pixel with an extreme 
intensity (maximum or minimum) in a region the size of the expected size of the 
small feature or object. The region should have the same number of dimensions as 
the data array, just fewer pixels. In other words, the region is defined as a sub- 
array of the multidimensional size equal to the expected size of the feature. In the 

25 case of mammograms, this sub-array is a square that is about 100 pixels in x and 
100 pixels in y when the resolution of the image is about 25 microns per pixel. To 
avoid selecting local extremes that are insignificant, the extreme is also required 
to achieve a certain absolute value - above a pre-set bright threshold in the case 
of a local maximum, or below a pre-set dark threshold in the case of a local 

30 minimum. 
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A pixel P on a ray is considered to be on the edge of a small object if it 
provides a maximum edge metric in a line search on a ray originating from the 
local extreme pixel and moving in a direction k. The edge metric may be defined 
as the change in intensity with each succeeding pixel in the direction k or by a 
5 Sobel operator centered on the pixel, or by any known edge metric. However, in 
the preferred embodiment with a local maximum, a ray edge pixel is found that 
more closely corresponds to that selected by expert analysis when the edge metric 
is a slope defined according to equation 3. 

,/», . f(p.)- f(p) 

S(P> - d(P„,P) P> 

10 For each pixel P around this local maximum P 0 the slope has a value S(P) where 
f(P 0 ) is the intensity, e.g., the gray scale value, at the local maximum pixel P 0 , and 
f(P) is the intensity at pixel P, and d(P 0 , P) is the distance between the local 
maximum pixel P 0 and the pixel P. In general, to extend to the case where P 0 is a 
local minimum, the absolute value of the numerator is used. The notation d(P„P 2 ) 

15 here indicates the absolute value of the distance between two points P, and P 2 . Let 
P n represent the nth pixel along a ray in a direction k. The n varies from 0 at the 
local maximum to N-l at the Nth consecutive pixel along the ray. The number N 
is not a critical choice as long as it is larger than the number of pixels expected to 
lie between the local maximum and the edge of the largest structures of interest. 

20 Referring to Figure 3, N should be the number of pixels extending half the length 
of the arrow 330 indicating the maximum expected size of a small feature, for 
example. Among the pixels P n , the pixel at which S(PJ is maximal is considered 
to be an edge point in that direction and is denoted by e(k). In the preferred 
embodiment, the ray search is applied in many equally spaced directions 

25 originating from the local maximum pixel, resulting in a set of ray edge pixels 
e(k) where k varies from 1 to K, the number of directions for which rays are 
computed. In the preferred embodiment, as shown in Figure 3, K equals 16. For 
each direction k, the edge pixel and all pixels between the local maximum and the 
edge pixel e(k) are labeled as belonging to the object associated with the local 
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maximum pixel P 0 . This results in K radial lines of labeled pixels 350, as shown 
in Figure 3. These labeled pixels are used as seeds or reference pixels for growing 
a region to identify all the pixels of the object. 

To identify all pixels lying within a contour including the edge points e(k), 
5 the region should grow essentially on pixels with more extreme intensity (e.g., 
increasing intensity) and toward the local extreme (e.g., local maximum). From 
any labeled pixel taken as a reference point, the region can grow to an adjacent 
unlabeled pixel if this new pixel satisfies some particular conditions. In the case 
of data arrays with more than two dimensions, adjacent points to a labeled point 

10 are those whose indices are all within one of the corresponding indices of the 
labeled point. Referring to Figure 4, the reference pixel is the labeled pixel 420 
and the eight adjacent pixels are numbered clockwise from the diagonally upper 
left pixel as pixel 1 through 8. These eight pixels are considered eight-connected 
with the labeled pixel 420. A subset of these adjacent pixels is the four-connected 

15 set of pixels to which pixels labeled 2, 4, 6 and 8 belong. With respect to the 
reference or labeled pixel 420, an eight-connected adjacent or neighbor pixel is 
checked. If the neighbor pixel is already labeled, it has already been determined 
that the neighbor pixel is on the object. If the neighbor pixel P is not labeled, then 
it has to satisfy the following conditions to be labeled. 

20 IF f(P)>f(P r ) then P must be in a position that constitutes a step from P r 

toward P 0 . 

IF f(P)<f(P r ), then P should be closer to P 0 than P r is to P 0 by more than a 
minimum distance called an inclusion tolerance distance. 

All pixels labeled during the process are used as reference pixels. The 
25 method stops when no pixel can be appended as shown in step 270 of Figure 2A. 
The step for labeling unlabeled pixels is illustrated in Figure 2A as step 260. 

The intensity and distance criterion referred to in step 260 are now 
described with reference to Figures 2C and 2D, which each show one of the two 
alternative criteria used in the present hill climbing method and apparatus. In both 
30 these figures, the first condition checked is the intensity f(P) of the unlabeled 
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point P compared to the intensity f(P r ) at the reference pixel P r , as shown in step 
262. 

Most microcalcifications have an intensity that decreases monotonically 
from the local maximum toward the edges. However, in some cases, this may not 
5 be true, and the growth toward the local maximum may need to include new 
pixels that have lower values or less extreme values than their labeled referenced 
pixels. As long as this is done strictly toward the local extreme, growth in an 
unwanted direction is avoided. That is, if the unlabeled pixel P is much closer to 
the local maximum (or minimum) than is the labeled referenced pixel P r , then the 

10 unlabeled pixel P is considered engulfed by the object and is labeled even if its 
intensity f(P) is less extreme than f(P r ). The distance by which the unlabeled point 
must be closer than the labeled point to be engulfed by the object is called the 
inclusion tolerance distance. In this and the following discussions, the difference 
in distances between the labeled and unlabeled points to the local maximum P 0 is 

1 5 represented by G given in Equation 4. 

G = d(P 0 ,P)-d(P 0 ,P r ) (4) 
When the unlabeled pixel P is closer to the local maximum P 0 than the unlabeled 
pixel ? 0 then G is negative. Therefore, the negative of G is compared to the 
inclusion tolerance to determine if the unlabeled pixel is close enough to the local 

20 extreme to be engulfed, as shown in step 263 of Figures 2C and 2D. In the 
preferred embodiment, the inclusion tolerance is one pixel. Thus, lower intensity 
pixels closer to the local maximum than the already labeled point P r by more than 
one pixel are close enough to be labeled. That is, a new pixel P with intensity f(P) 
less extreme than the intensity f(P r ) of the referenced pixel P r is appended to the 

25 region if its distance to the local extreme is such that -G is > the inclusion 
tolerance distance, as shown in step 265 of Figures 2C and 2D. If the unlabeled 
pixel with less extreme value is less than the inclusion tolerance closer to the local 
extreme or is farther from the local extreme, then the unlabeled pixel is not 
labeled, as shown in step 265 of Figures 2C and 2D. 
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The other branch from step 262 in Figures 2C and 2D is followed when 
the adjacent pixel P that is unlabeled has an intensity that is greater than or equal 
to the intensity of the labeled pixel P r This corresponds to the condition in the 
case of a local minimum that the unlabeled pixel has a lower intensity than the 
5 labeled pixel P r . That is, the << y es " branch is followed from box 267, in general, if 
the unlabeled pixel P has an intensity that is no less extreme than the intensity at 
the labeled pixel P r . Each of two different criteria can be used to determine 
whether the unlabeled pixel P is in a position that constitutes a step from the 
labeled pixel P r toward the extreme pixel P 0 . 

10 The first criterion, Criterion 1, is indicated in Figure 2C and step 264a and 

is based on the angle of the line perpendicular to the line segment connecting the 
local extreme P 0 with the reference pixel P r The line perpendicular to the segment 
connecting the local extreme to the labeled pixel is called the reference line 430 
and is shown in Figure 4. For arrays of more than two dimensions, the reference 

15 would be a surface with a number of dimensions at least one dimension less than 
the multidimensional array. The numbered pixels of Figure 4 are approved for 
appending to the small feature if they fall within the list of approved pixels listed 
in Table 1 for the quadrant in which the angle 9 varies from 0-90°. The angle 9 
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between the reference line 430 and the Jt-axis is also shown in Figure 4. The first 
two columns of Table 1 show the relationship between the coordinates of the 
reference pixel x, and y r of P r and their relationship to the coordinates Xq and y 0 of 
the local maximum P 0 . For different values of the angle 0 or its tangent, tan 9, 
5 different of the numbered pixels in Figure 4 are approved. Table 1 captures the 
condition that the unlabeled pixel P and the local maximum P 0 must lie on the 
same side of the reference line 430. Among the eight pixels that surround a 
reference pixel, only some will meet the spatial criterion of Criterion 1, depending 
on the angel G of the reference line. The angle G is measured positive 
10 counterclockwise from the *-axis. The allowable pixels for values of G in the other 
three quadrants are obtained in a symmetrical manner. An extended table would 
have to be drafted for data arrays of greater than two dimensions. 

As an alternative for the constraint CI described above and summarized in 

15 Table 1, Constraint 2 can be used to determine whether a neighboring pixel 
should be labeled. Constraint 2 is more readily extensible to more than two 
dimensions. Referring to Equation 4 defining the distance difference G, most 
allowable pixels described by Criterion 1 yield a negative G value. However, 
some pixels generate a positive G value. These positive G pixels are the pixels 

20 that provide a step, from the reference pixel P r , approximately parallel to the 
reference line. This type of growth through pixels is especially desirable around 
the edge of the small structure. The largest values of G are associated with 
diagonal pixels and occur at the edge of the smallest features to be segmented. 
Furthermore, among all possible pixel configurations, the value of G is maximal 

25 when the reference line angle Q is 45° or 135° and the new pixel P is diagonally 
connected to the referenced pixel P r . This maximal value is also obtained for other 
homologous arrangements of the three pixels. A positive threshold G t for G can be 
used instead of Criterion 1. Consider an approximately circular object 2N pixels 
wide. On the edge of such an object, the highest value for G, called G^, will 

30 equal (V(N 2 +2))-N. The smaller N, the larger G^ will be. An appropriate 
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threshold for G can be set by using the width of the smallest object of interest. 
Therefore, an alternative way of constraining the expansion of pixels away from 
the local extreme is to allow only new pixels which provide a value of G of at 
most G^. That is, set G^G^. This threshold, G t , can be considered an expansive 
5 tolerance distance. Criterion 2 can be stated as: G must be less than or equal to the 
expansive tolerance distance, G t . For example, mammograms with pixels of 25 
microns and microcalcification candidates having structures as small as 0.25 mm 
across, yield N=5; so, G t = G^ = 0.196. 

The preferred embodiment determines 16 ray edge pixels around the 

10 object, and segments with the hill climbing procedure described. As indicated in 
step 270 of Figure 2A, each appended pixel is labeled and is used as a reference 
pixel itself during growth. The growth stops when no pixel can be appended. 
Once no more new pixels can be labeled, each labeled pixel is examined to 
identify edge pixels of the small feature in step 275 in Figure 2B. The edge pixels 

15 of the small feature are determined to be all labeled pixels that are four-connected 
to an unlabeled pixel after no further pixels can be added. 

After every object has been segmented and its outer edge pixels defined, 
larger features may be discernable. The larger features can be constructed where 
the small features abut or overlap slightly. The step of joining small features 

20 together into a larger feature is depicted in step 285 of Figure 2B. Depending on 
the larger feature being assembled, the criterion for joining small features can be 
that the small features share edge pixels, or that the edges overlap so that the edge 
of one small feature is an interior labeled pixel of another small feature. It is also 
possible that features be joined that do not touch or overlap, provided they are 

25 sufficiently close together. A tolerance called a join distance can be used to 
determine how close the edges should be to each other in order to combine the 
small features into one or more larger features. In this case, all small features are 
joined where the edge pixels of two different small features are within a join 
distance. Overlapping pixels are covered by this criterion as are features whose 
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edge pixels coincide. By setting the joined distance to 0 the edge coincidence is 
required; and by setting the join distance negative, overlapping can be required. 



EXAMPLES 

To determine whether the results of the present invention provide edges of 
5 small features that are useful in interpreting mammograms and in doing so with 
fewer computations than other methods, several experiments were performed with 
actual mammograms. The correctness of the edge determined by the present 
invention is measured by its similarity to the edges determined by an analyst, and 
its ability to discriminate among the candidate microcalcifications in subsequent 
10 processing. Other advantages of the preferred embodiment are measured using the 
complexity or number of computations involved in the procedure, and the time 
required to execute the procedure on a computer. 

Example 1 

Five mammograms containing subtle microcalcification clusters were used 
15 to evaluate the algorithms for data that would warrant the use of an automated 
system. Mammograms without magnification were used; and the breast images 
covered an area that ranged between 12cm x 6cm and 21cm x 1 1cm. The location 
of individual microcalcifications was indicated by an experienced 
mammographer. These 5 mammograms contained 15 clusters with a total of 124 
20 microcalcifications, yielding about 8 microcalcification per cluster. The number 
of microcalcifications per cluster ranged between 3 and 18. The size of 
microcalcifications ranged between 0.25mm and 1mm wide, with more than 90% 
being smaller than 0.5mm. Mammograms were digitized with a Howtek D4000 
drum scanner using a spatial resolution of 25 microns per pixel and 12-bit A/D 
25 conversion, with an optical dynamic range of 0-3.5 optical depths (O.D.). 

The multi-tolerance region growing procedure grows a region around a 
seed pixel by appending 4-connected pixels P that satisfy: 

(l^)(F_+F m J/2>P>(l-x)(F rmx +F min )/2 (5) 
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where t is the tolerance parameter, and and F min are the current maximum and 
minimum values in the region grown that far. The value of x is not manually 
selected by the user; the best value is automatically determined for each 
segmented structure by repeating the growth with multiple values of t between 
5 0.01 and 0.4 with steps of s=l/v, where v is the 8-bit value of the seed pixel. 
Three features are extracted from each region grown with a different tolerance 
level: shape compactness, center of gravity, and size. The algorithm determines 
the value of x that results in the minimal change in the vector of these three 
features with respect to the previous x value in the sequence by computing a 

10 normalized difference between consecutive vectors. The vector with minimal 
difference indicates the best choice of x. 

The segmentation outcome of the multi-tolerance region growing 
procedure on 5 subtle microcalcification candidates depended partly on the 
intensity structure of the microcalcification. When the intensity transition from the 

15 edge to the background was relatively abrupt, the segmented region coincided 
closely to the visually perceived edge. When the intensity at the edge decreased 
gradually toward the background level, this algorithm generally produced a 
relatively large region. Nevertheless, the growth was consistently contained, i.e. it 
did not grow to an unacceptable size and it generated boundaries that can be used 

20 as an estimate of the immediate background around the microcalcification. 

The active contours model represents the contour points as v(s)=(x(s),y(s)) 
The contour is obtained by minimizing the energy functional: 



where E im is the internal energy due to the elasticity and the rigidity, PE is the 
25 potential energy obtained from the image data, E cxt is the energy of external forces 
that can be applied to the contour. The integration is performed over the entire 
contour Q. The internal energy is expressed by: 




(6) 



E int =w 1 |v'(s)| 2 +w 2 |v"(s)| 2 



(7) 
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where w, and w 2 are coefficients that control the elasticity and rigidity, 
respectively, and primes denote differentiation. The choice of potential energy 
depends on the application; it is typically the negative squared gradient 
magnitude, and is so used for mammograms. 
5 The active contour that minimizes E(v) satisfies the Euler-Lagrange 

equation: 

-(w I v')'+(w 2 v'T=F(v) (8) 
where F(v) represents the force due to the combined effects of the potential energy 
and external energy. In this study we implemented the balloon forces and the 
10 image force normalization suggested, resulting in 

F(v)=k in (s)-k 2 ™| (9) 



where n(s) is the unit vector normal to the contour at point v(s), oriented toward 
the outside of the contour, k, is the magnitude of the balloon inflation force, and 
k 2 is the coefficient of the normalized image force. The value of k 2 is selected to 

15 be slightly larger than k, to allow edge points to stop the inflation force. 

The numerical solution was implemented using finite differences and the 
iterative evolution as suggested: 

a+rA)v=(v M -HuF(v M )) (10) 
where I is the identity matrix, x is the time step, A is the pentadiagonal matrix 

20 obtained with the finite difference formulation of E int , v t is the active contour 
vector at time t, and F(v t ) is the external force vector at time t. We used the 
negative squared magnitude of the image gradient as the potential energy. Pixels 
detected with an edge detector were not used in this study. The gradient of the 
image was computed with the Sobel operator. 

25 The initial position of the contour was set automatically for each structure 

to be segmented. Since each structure of interest is a local intensity extreme, 
pixels were selected that were local maxima across the entire image. Each local 
maximum was used to segment a region around it. The width of the smallest 
microcalcifications considered in this study was about 0.25 mm and the majority 
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of the microcalcifications in our database had widths in the range 0.3 to 0.5 mm. 
A circle of 0.2 mm diameter around the local maximum pixel was used as the 
initial position of the active contour. The initial contour points were 24 8- 
connected pixels forming this circle. 
5 The selection of parameters for the active contour segmentation required 

some trial and error to obtain good segmentation. The segmentation of the same 5 
subtle microcalcification candidates was performed using different active contours 
parameters. First, following the recommendations of Cohen (Cohen, L.D. "On 
Active Contour Models and Balloons," CV GIP: Image Understanding, vol. 53, 

10 pp. 211-218, 1991), we selected the values of w, and w 2 as a function of the 
spatial discretization step size h, such that w, was of the order of h 2 and w 2 was of 
the order of h 4 (w,=6, w 2 =40). Then x was also set to 0.1. When k x and k 2 were 
relatively small (2 and 4), the image force and the balloon force did not act 
sufficiently on the active contour, producing contours that were only slightly 

15 different than the initial position. When these two parameters were increased (14 
and 16), the resulting segmentation was very close to that expected visually. 
Increasing these parameters further (24 and 26) increased the combined effect of 
image gradient and balloon forces, producing contours that extended beyond the 
expected edges. Within this range, segmentation with the active contour model 

20 was not very sensitive to the values of the other parameters. The effect of 
doubling w, to 12, is that contours became slightly smaller due to the increased 
stiffness of the active contour model. Sensitivity to w 2 was also low. When w 2 was 
doubled to 80, the contours became slightly smoother due to the increased rigidity 
of the model. 

25 The segmentation steps of the hill climbing approach of the present 

invention are illustrated in Figure 6. Figure 6A shows a microcalcification 
candidate that has a width of about 0.3 mm. The 16 ray edge points 624 
determined by the radial line search of the hill climbing algorithm are shown in 
Figure 6B. The region grown using spatial Constraint 1 is in Figure 6C. The 

30 region grown with spatial Constraint 2 was identical for this microcalcification 
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candidate. The edge pixels 642 of the entire microcalcification candidate are 
shown in Figure 6D. The segmentation of microcalcifications by the hill climbing 
method produces outcomes using the spatial Constraints 1 and 2 that were almost 
identical. In this study, about a quarter of microcalcifications were segmented 
5 identically by the two spatial constraints and the rest differed by a few pixels, 
resulting in a negligible change over the entire microcalcification. Both spatial 
constraints directed the growth of the regions successfully, resulting in regions 
that were compatible with visual interpretation. 

The differences between the three methods are illustrated in Figure 7. 

10 Three subtle microcalcifications candidates are shown in Figure 7A. When the 
contrast of a microcalcification candidate was relatively low, or parts of it 
exhibited a very gradual decrease in intensity toward the background, the multi- 
tolerance algorithm (Figure 7B) segmented a larger region than those of the other 
two algorithms. Good segmentation with active contours (Figure 7C) was 

15 obtained using Wj=6, w 2 =40, t=0.1, k,=14 and k 2 =16, for all microcalcifications 
candidates of this study. Using these parameters, segmentation with active 
contours provided edges 735 that were smoother than edges 725 and 745 
produced by segmentation with the other two methods. The selection of Wj and w 2 
provided the flexibility needed to adapt relatively well to the shape of diverse 

20 microcalcifications candidates. The elasticity level allowed the contour to grow to 
the highest gradient locations when the segmented structures were relatively large, 
and the rigidity level allowed the contour to develop sharp bends dictated by the 
data in some microcalcifications. The edges 745 of regions grown by the hill 
climbing algorithm shown in Figure 7D were not as smooth as those 735 of the 

25 active contours, but the convolutions were consistent with visually perceived 
edges around microcalcifications candidates. 



Example 2 

Segmentation of microcalcification candidates serves as an initial step for 
discriminating between the population of microcalcifications and that of 
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background structures. The discrimination potential of each segmentation 
algorithm was quantified using features extracted from structures segmented 
around all the local maxima in the 5 mammograms. These structures consisted of 
the 124 microcalcifications mentioned above and 2,212 background structures 
5 segmented in the same mammograms. Four characteristics were used to assess the 
discrimination potential in this study. 

1. Contrast was measured as the gray level difference between the 
local maximum pixel P 0 in the structure, and the mean of pixels around its edge. 

2. Relative contrast was computed as the ratio of the contrast to the 
10 value at the local maximum. 

3. Area was computed as the number of labeled pixels in the grown 

region. 

4. Edge sharpness was the mean of the gradient computed with a 
Sobel operator across all edge pixels. The Sobel operator is a mask which weights 

15 the eight neighbors of a pixel to compute a sum proportional to the gradient x, or 
the y gradient, or total gradient. 

The discrimination ability of each feature was determined separately using 
the area under a receiver operating characteristic (ROC) curve obtained with that 
feature. The ROC curve pots the percentage of correctly detected 

20 microcalcifications against the percentage of detected background structures as a 
detection threshold is changed. The ROC curve area is higher when the feature 
has distributions that are more separable for a given property. When both 
populations overlap completely, the ROC curve area is 0.5. In general, effective 
discrimination power is indicated by a value above 0.8. Table 2 summarizes the 

25 results for all three procedures. The area feature had very low discrimination 
power for all three algorithms, indicating that the two types of structures cannot 
be discriminated well on the basis of their area segmented. However, the other 

Table 2 

Multi-tolerance Active Contours Hill Climbing 
Region Growing 
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Contrast 



0.80 



0.82 



0.83 



Relative Contrast 



0.83 



0.90 



0.90 



Area 



0.63 



0.60 



0.54 



Sharpness 



0.80 



0.85 



0.85 



three features suggested good discrimination potential for all three algorithms. A 
comparison among algorithms shows that both the hill climbing method of the 
present invention and the active contours algorithm provide segmentation with the 
same discrimination power, and they both perform slightly better than the multi- 
tolerance segmentation. Thus, the hill climbing method produces edges as good as 
the best produced by the conventional approaches tested. 

The significant advantage of the hill climbing algorithm is its speed. While 
the multi-tolerance algorithm provides a good solution to avoid the use of 
statistical models, local statistics estimators and the manual selection of threshold, 
its cost is multiple segmentations of the same structure and computation of 
features during the segmentation of each structure. Furthermore, in some cases, 
this algorithm segments regions that are somewhat larger than expected. 
Consequently, the time required for segmentation of a mammogram with this 
algorithm is high. The segmented regions were comparable to those of the other 
two algorithms in many cases. The differences were caused by the fact that the 
growth mechanism of this algorithm is constrained only by an intensity range 
criterion applied to a new pixel. In contrast, active contours are constrained by 
internal forces that regulate the growth away from the local maximum, and hill 
climbing has an inward growth mechanism based on edge points. 

The active contours also circumvent the statistical and manual threshold 
selection issues for each mammogram, but the selection of the operational 
parameters for a set of mammograms requires some trial and error. However, 
when an appropriate set of parameters is determined, it appears to be valid for a 
wide range of microcalcifications so it need not be modified with each 
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mammogram. The choice of negative squared gradient magnitude as the image 
energy function seems to be an appropriate one to segment microcalcifications. 

Example 3 

The computational complexity c ra of the multi-tolerance region growing 
5 algorithm is of the order 0(4smo) where s is the number of steps in the tolerance 
search, m is the number of pixels in the region, and o is the number of operations 
per pixel. The factor 4 is included because the algorithm visits the 4-connected 
neighbors for each pixel in the region. Considering 125 to be an average intensity 
value for the local maximum, the average step size is 0.008 resulting on the 

10 average in about s=50 steps to cover the range 0.01 to 0.4. The average size of 
segmented structures is about 200 pixels. At each pixel the computations 
performed include intensity comparisons, update of and F^, and calculation 
of the center of gravity. Considering about 12 operations per pixel on the average, 
the numerical estimate for the average number of operations per segmentations is 

15 c m =480,000. 

The computational complexity c a of the active contour model is 
0[2(n+n 2 )t] where n is the number of contour points, and t is the number of 
iterations. The factor of 2 is included due to the fact that the x and y coordinates 
of each contour point are computed separately, with identical operations. At each 

20 iteration, order n computations are needed to determine the normal vectors, and 
order 2n 2 operations are needed to perform a matrix multiplication. In this study 
24 contour points were used, and the number of iterations depended on the size of 
the structure. On the average however, the active contour model converged in 
about 20 iterations. This resulted in an average value of c a =47,040, a factor of ten 

25 improvement over the multi-tolerance method. 

The complexity c h of the hill climbing method is 0(KN +8m) where K is 
the number of radial directions from the local maximum, N is the number of 
pixels searched in each direction, and m is the number of pixels in the grown 
region. A factor of 8 is included since all 8 neighbors of each pixel are visited. In 
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this study K was 16 and N was 40, and considering an average structure size of 
m=200, the average estimate of the number of operations is (^=2,240, a factor of 
20 improvement over the active contour methods, and 200 over the miilti- 
tolerance method. The proportions of c m , c a and are approximately 214:21:1 
5 respectively, with hill climbing far less complex than the other two methods. 

Example 4 

The speed of the different methods was compared using a section of a 
mammogram containing 456 local maxima, 35 of which were in 
microcalcifications. The sizes of microcalcifications ranged between 0.25 mm and 

10 0.5 mm. The times to complete the segmentation of this section of mammogram 
using the three algorithms implemented in C on a 10 million floating point 
operations per second (MFLOPS), IBM 6000 computer were 17 minutes 47 
seconds for the multi-tolerance algorithm, 1 minute 47 seconds for the active 
contours, 7 seconds for hill climbing with spatial Constraint 1, and 5.4 seconds for 

1 5 hill climbing with spatial constraint 2. 

Hill climbing with spatial Constraints 1 and 2 yielded practically identical 
segmentations; but the method was about 20% faster using spatial constraint 2, 
resulting in 1 1.8 ms on average for segmenting a structure, as opposed to 15.3 ms 
obtained with spatial Constraint 1. 

20 A common technique to determine the edges of an object uses an edge 

enhancement algorithm such as the Sobel operator, thresholding to separate the 
pixels on edges, and pixel linking to string edge pixels that belong to the same 
object. Selection of the threshold is critical, and linking poses problems in 
segmenting microcalcifications because there are many closely spaced small 

25 structures in a background that are likely to produce considerable numbers of edge 
pixels. The hill climbing method of the preferred embodiment determines edge 
points that are on the edge of the same object by virtue of the radial line search 
emanating from the same local maximum. It does not require a threshold to 
separate edge pixels because the slope in Equation 3 is referred to the local 
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maximum and is greatest at pixels that are on, or very near, the visually perceived 
edges. Finally, the hill climbing method avoids some pitfalls of the region 
growing mechanism by growing a region inward, toward the local maximum. 

There has been disclosed a segmentation method and apparatus for data 
arranged in a multidimensional array which overcomes the problems of the prior 
art. Although the present invention has been described above by way of detailed 
embodiments thereof, it is clearly understood that variations and modifications 
may be made by one of ordinary skill in the art and still lie within the spirit and 
scope of the invention as defined by the appended claims and their equivalents. 



